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Abstract 

Biological functions in living cells are controlled by protein interaction and genetic 
networks. These molecular networks should be dynamically stable against various 
fluctuations which are inevitable in the living world. In this paper, we propose and 
study a stochastic model for the network regulating the cell cycle of the budding 
yeast. The stochasticity in the model is controlled by a temperature-like parameter 
(3. Our simulation results show that both the biological stationary state and the 
biological pathway are stable for a wide range of "temperature". There is, how- 
ever, a sharp transition-like behavior at /3 C , below which the dynamics is dominated 
by noise. We also define a pseudo energy landscape for the system in which the 
biological pathway can be seen as a deep valley. 
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1 Introduction 

Quantitative understanding of biological systems and functions from their 
components and interactions presents a challenge as well as an opportunity 
for interested scientists of various disciplines. Recently, a considerable amount 
of attention has been paid to the quantitative modeling and understanding of 
the budding yeast cell cycle regulation [1,2,3,4,5,6,7,8,9,10]. In particular, Li et 
al. [3] introduced a deterministic Boolean network model and investigated its 
dynamic and structural properties. Their main results are that the network is 
both dynamically and structurally stable. The biological stationary state is a 
global attractor of the dynamics; the biological pathway is a globally attracting 
dynamic trajectory. These properties are largely preserved with respect to 
small structural perturbations to the network, e.g. adding or deleting links. 
However, one crucial point left unaddressed in their study is the effect of 
stochasticity or noise, which inevitably exists in a cell and may play important 
roles [11]. In this paper, we advance a probabilistic Boolean network [12,13] 
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Fig. 1. The cell-cycle network of the budding yeast. Each node represents a protein or 
a protein complex. Arrows are positive regulation, "T" -lines are negative regulation, 
dotted "T" -loops are degradation. 

on the protein interaction network of yeast cell cycle. We found that both 
the biological stationary state and the biological pathway are well preserved 
under a wide range of noise level. When the noise is larger than a value of 
the order of the interaction strength, the network dynamics quickly becomes 
noise dominating and loses its biological meaning. 



2 Method 



Our stochastic model is based on the updated protein interaction network of 
Li et al. [3], which is depicted in Fig. 1. Nodes in the figure represent proteins 
or protein complexes. Arrows represent positive interaction, or "activation". 
Lines with a bar at the end represent negative interaction, or "repression". 
Dotted loops with a bar represent self-degradation. We refer the reader to 
Ref. [3] for a full biological account of this network. Here we only give a very 
brief summary. There are four phases in the cell cycle process: the Gl phase 
in which the cell grows, the S phase in which the DNA is copied, the G2 
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phase in which the cell prepares for mitosis, and the M phase in which the 
two chromosome copies are separated and the cell divides into two. There are 
several checkpoints during the process to ensure that the next event will not 
happen until the current event is finished. So the process could be blocked at 
checkpoints. Following Ref. [3], we keep only one such checkpoint in the model: 
"cell size" . Thus the picture for the cell division process is the following: The 
cell is resting on a stationary state Gl (blocked at the checkpoint until it 
grows big enough). The "signal" to start the cell-cycle process comes from the 
"cell size" which turns on a cyclin Cln3. Cln3 activates a pair of nodes, SBF 
and MBF. SBF and MBF stimulate the transcription of Gl/S genes, including 
those of Cln2 and Clb5. The S phase cyclin Clb5 initiates DNA replication, 
after which the transcription factor complex Mcml/SFF is turned on, which 
stimulates the transcription of many G2/M genes including the gene of the 
mitotic cyclin Clb2. The cell will exit from mitosis and divide into two after 
Clb2 is being inhibited and degraded by Cdc20, Cdhl and Sicl. The cell (or 
two cells: the mother and the daughter) now comes back to the stationary Gl 
state, waiting for the signal for another round of division. So from a dynamics 
point of view, the cell's stationary state Gl is a fixed point. A "start" signal 
will take it out of the fixed point, and it will then go through a specific dynamic 
trajectory (the biological pathway for cell division), and come back to the fixed 
point. 

In our model, the 11 nodes in the network shown in Fig. 1, namely, Cln3, MBF, 
SBF, Cln2, Cdhl, Swi5, Cdc20, Clb5, Sicl, Clb2, and Mcml are represented 
by variables (s±, s 2 , Sn), respectively. Each node % has only two values: 
Si = 1 and Si = 0, representing the active state and the inactive state of the 
protein i, respectively. Mathematically we consider the network evolving on 
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the configuration space S = {0, 1} ; the 2 11 = 2048 "cell states" are labelled 
by {n = 0, 1, ...,2047}. The statistical behavior of the cell state at the next 
time step is determined by the cell state at the present time step. That is, the 
evolution of the network has the Markov property [14] . The time steps here are 
logic steps that represent causality rather than actual times. The stochastic 
process is assumed to be time homogeneous. Under these assumptions and 
considerations, we define the transition probability of the Markov chain as 
follows: 

P r ( Sl (t + 1), s u (t + l)|si(f), Sn(f)) 
n 

= n^r(Si(*+l)|si(*),...,Sn(*)), (1) 

i=l 

where 

P r (s i {t + l)=a i \s 1 (t),...,s 11 (t)) = 

exp(/3(2a, - 1)T) 
exp(/3T) +exp(-/3T)' 

if T = E)Li a ijSj (t) ^O.a.e {0, 1}; and 

P r ( Si (t + 1) = Si (t)\ Sl (t), s n (t)) = T-^, (2) 

1 + e a 

if T = Y}jL\C'ijSj{t) = 0. We define = 1 for a positive regulation of j 
to i and a^- = —1 for a negative regulation of j to i. If the protein i has a 
self-degradation loop, an = —0.1. The positive number j3 is a temperature- 
like parameter characterizing the noise in the system [15]. Noticeably, the 
actual noises within a cell might not be constant everywhere, but here we use 
a system-wide noise measure for simplicity. To characterize the stochasticity 
when the input to a node is zero, we have to introduce another parameter a. 
This parameter controls the likelihood for a protein to maintain its state when 
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there is no input to it. Notice that when (3, a — > oo, this model recovers the 
deterministic model of Li et al. [3]. In this case, they showed that the Gl state 
(the purple node in Fig. 3) is a big attractor, and the path (blue nodes — ► 
olive-green nodes — > dandelion nodes — > red nodes — > purple node in Fig. 3) is 
a globally attracting trajectory. Our study focuses on the stochastic properties 
of the system. 

Because the Markov chain consists of finite states and is irreducible, every 
state is accessible to all others. Therefore all of the 2048 states constitute 
a communicated recurrent class, and the Markov chain is ergodic. In this 
case, there exists a probability distribution II = (7r , 7Ti, 7^047), an invariant 
measure, such that for all states m,n e {0, 1, ...,2047}, 

lim p mn {r) = TT n 

where p mn (r) is r-step transition probability from the initial state m to the 
target state n. That is to say, when r is big enough, the probability for the 
system to reach state n is almost independent of the starting position m. Even 
though each state has a positive probability, the order of the magnitudes of 
the probabilities are very different among the states; some are so small that 
in realistic cases they can never be observed. 

Our interests are in the asymptotic behavior of the dynamic system. Steady 
state probability distribution II = (tto, tti, 7^047) can be found by solving 
linear equations IIP = II, where P is the transition matrix of the Markov 
chain. The net probability flux from m to n is then n m p mn — TT n Pnm, where 
p mn is the transition probability from m to n. For a given (3, one can define a 
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pseudo energy for the state n as [16] 



E n — 



log7r 5 



n 



(3) 



We first study the property of the biological stationary state Gl and define an 
"order parameter" as the probability for the system to be in the Gl state, ttqi- 
Plotted in Fig. 2 is the value of the order parameter (ttgi) as a function of the 
control parameter (3 with different a. At large (3 (low "temperature" or small 
noise level), the Gl state is the most probable state of the system and 7r G1 has 
a significant value. Note that for a finite a there are always "leaks" from the 
Gl state, so that the concept of attractor in the deterministic model in Ref. 
[3] cannot applied here, hqi decreases with the decrease of (3; one observes a 
transition-like behavior as the function of (3 (similar behavior has been seen 
in [13]). In order to compare this transitional behavior to the transition in the 
system of thermodynamic equilibrium, we define 



to fit the order parameter (ttgi) curves in Fig. 2, where (3 C , b, a are parameters. 
When a is fixed to 5, we obtain (3 C « 1.03, b « 0.36, a « 0.36 (See Fig. 2(b)). 

At around (3 C = 1.03, hq\ drops to a very small value, indicating a "high 
temperature" phase in which the network dynamics cannot converge to the 
biological steady state Gl. The system is, however, rather resistant to noise. 
The "transition temperature" is quite high: the value of (3 C ~ 1.03 implies that 
the system will not be significantly affected by noise until approx 10% of the 
updating rules are wrong (e~ im / (e im + e~ L03 ) 0.1). 

We next study the statistical properties of the biological pathway of the cell- 



m = b\ 




a 



(4) 
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Fig. 2. The probability of the stationary Gl state and the biological pathway. The 
order parameter ttgi as a function of f3 with (a), a = 2; (b), a = 5, and (c),a = 10. 
The solid line in (b) is the fitting function tp(/3) = 0.36[y^3 — 1| - 36 . 




Fig. 3. The probability flux. For a node, only the largest flux from it is shown. 
The nodes on the biological pathway are denoted with different colors: purple, the 
stationary Gl state; blue, the other Gl states; olive-green, S state; dandelion, G2 
state; and red, M states. All other states are noted in normal green. The simulations 
were done with a = 5 and (3 = 6. 

cycle network. We search the probability for the system to be in any of the 
biological states along the biological pathway, as a function of (3. One observes 
a similar transition-like behavior as shown in Fig. 2. The jump of the proba- 
bility of the states along the biological pathway in the low temperature phase 



8 



is due to the fact that in this phase the probability flux among different states 
in the system is predominantly the flow along the biological pathway. To visu- 
alize this, we show in Fig. 3 an example of the probability flux among all 2048 
states. Each node in Fig. 3 represents one of the 2048 states. The size of a node 
reflects the stationary distribution probability of the state. If the stationary 
probability of a state is larger than a given threshold value, the size of the 
node is in proportion to the logarithm of the probability. Otherwise, the node 
is plotted with the same smallest size. The arrows reflect the net probability 
flux (only the largest flux from any node is shown). The probability flux is 
divided into seven grades, which are expressed by seven colors: light-green, ca- 
nary, goldenrod, dandelion, apricot, peach and orange. The warmer the color 
is, and the wider the arrow is, the larger the probability flux. The width of 
an arrow is in proportion to the logarithm of the probability flux it carries. 
The arrow representing the probability flux from the stationary Gl state to 
the excited Gl state (the START of the cell-cycle) is shown in dashed lines. 
One observes that once the system is "excited" to the START of the cell cycle 
process (here by noise (a) and in reality mainly by signals like "cell size") the 
system will essentially go through the biological pathway and come back to 
the Gl state. Another feature of Fig. 3 is that the probability flux from any 
states other than those on the biological pathway is convergent onto the bio- 
logical pathway. Notice that this diagram also characterizes the properties of 
fixed points that are ignored by Li (Ref. [3]). Those fixed points also converge 
onto biological pathway. For (3 < (3 C) this feature of a convergent high flux 
bio-pathway disappears. 

In the previous discussions, we see that there is a "phase transition" as a 
function of the "temperature" in the stochastic cell-cycle model. The next step 
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is to try to understand this transition-like behavior. For this purpose, we define 
a "potential" function and study the change of the "potential landscape" as 
a function of (3. Specifically, we define 

S n = - log vr„ = (3E n , (5) 



where E n is the pseudo energy defined in Eq. [3]. Fig. 4 shows four examples 
of AS n = S n — S0 distribution, where the reference potential SO in each plot is 
set as the highest potential point in the system. Note that the 11-dimensional 
phase space is reduced to two dimensions and there is no distance metric 
among the states in the 2D phase space. The states in 2D are arranged for 
clarity, which reflect a kind of dynamic relationship as in Fig. 3. One observes 
that far from the critical point ((3 = 0.01, Fig. 4(a)), the potential values 
are high (around -4), and the landscape is flat. Near but below the critical 
point (f3 = 0.6, Fig. 4(b)), some local minima (blue points) become more 
pronounced, but the landscape still remains rather flat. We notice that these 
minimum points do not correspond to the biological pathway. Right after 
the critical point (j3 = 1.5, Fig. 4(c)), the system quickly condenses into a 
landscape with deep valleys [17]. The state with the lowermost potential value 
corresponds to the stationary Gl state. A linear line of blue dots from up-left 
to down-middle corresponds to the biological pathway, which forms a deep 
valley. Some deep blue dots out of the biological pathway are local attractors 
in Ref. [3]. Notice that although their potential values are low, they attract 
only a few nearby initial states-all these points are more or less isolated. After 
the critical point, the potential landscape does not change qualitatively (see 
Fig. 4(d) with (3 = 6). As (3, a — > oo, the landscape becomes nine deep holes, 
each corresponding to an attractor of the determinate system. 
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Fig. 4. The "potential" landscape of the system before and after the critical point, 
(a) = 0.01, (b) = 0.6, (c) = 1.5, (d) = 6.0, all for a = 5. The color code 
gives the relative value of the potential function. 

3 Conclusion 

In conclusion, we introduced a stochastic model for the yeast cell cycle net- 
work. We found that there exists a transition-like behavior as the noise level 
is varied. With large noise, the network behaves randomly; it cannot carry 
out the ordered biological function. When the noise level drops below a criti- 
cal value, which is of the same order as the interaction strength (0 C « 1.03), 
the system becomes ordered: the biological pathway of the cell cycle process 
becomes the most probable pathway of the system and the probability of devi- 
ating from this pathway is very small. So in addition to the dynamical and the 
structural stability [3], this network is also stable against stochastic fluctua- 
tions. We used a pseudo potential function to describe the dynamic landscape 
of the system. In this language, the biological pathway can be viewed as a 
valley in the landscape [18,19,17]. This analogy to equilibrium systems may 
not be generalizable, but it would be interesting to see if one can find more 
examples in other biological networks, which are very special dynamical sys- 
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terns. 
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